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Abstract 

We study a Monte Carlo algorithm for simulation of probability dis- 
tributions beised on stochastic step functions, and compare to the tradi- 
tional Metropolis/Hastings method. Unlike the latter, the step function 
algorithm can produce an uncorrelated Markov chain. We apply this 
method to the simulation of Levy processes, for which simulation of un- 
correlated jumps arc essential. 

We perform numerical tests consisting of simulation from probability 
distributions, as well as simulation of Levy process paths. The Levy 
processes include a jump-diffusion with a Gaussian Levy measure, as well 
as jump-diffusion approximations of the infinite activity NIG and CGMY 
processes. 

To increase efficiency of the step function method, and to decrease 
correiations in the Metropolis/Hastings method, we introduce adaptive 
hybrid algorithms which employ uncorrelated draws from an adaptive dis- 
crete distribution defined on a space of subdivisions of the Levy measure 
space. 

The nonzero correlations in Metropolis/Hastings simulations result in 
heavy tails for the Levy process distribution at any fixed time. This prob- 
lem is eliminated in the step function approach. In each case of the Gaus- 
sian, NIG and CGMY processes, we compare the distribution at t = 1 
with exact results and note the superiority of the step function approach. 
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1 Introduction 



Levy processes are a type of stochastic process whose paths can behave erratic- 
ally like a Brownian motion, as well as include discontinuities, i.e. jumps. Many 
observed phenomena in nature and human society display continuous erratic mo- 
tion similar to Brownian motion, while also having random jumps at random 
times, and can therefore be modelled using Levy processes. Areas of applica- 
tion include e.g. microscopic physics, chemistry, biology and financial markets. 
Therefore the study of computer simulation methods for Levy processes is an 
important subject. 

In our case, will will employ a jump-diffusion approximation when we apply 
our method to infinite activity Levy processes. Jump-diffusion processes can be 
considered as a sum of three independent component processes. A deterministic 
drift, a Brownian motion, and lastly a finite activity jump process. The jump 
process is completely described by a Levy measure v on R. The average jump 
rate, or intensity, is given by the total weight A := and the distribution 

of jump values follows the Levy probability measure := vjX. 

By definition, the Brownian motion has independent increments, and this 
is also the case for subsequent jumps in the jump process. In computer sim- 
ulations, violation of these properties will give incorrect results. Simulation of 
the Brownian motion is easy, since well known algorithms exist for producing 
uncorrelated draws from a normal distribution, using a good pseudo-random 
number generator (PRNG). 

On the other hand, the simulation of the jump process requires more care. In 
principle we can obtain uncorrelated draws from v\ by inverting its cumulative 
distribution function. However, this might not be feasible to do for a given v\. 
In this case, the Metropolis/Hastings (MH) algorithm might come to the rescue. 
It can easily produce a Markov chain with values distributed according to v\. 
However, subsequent values will be correlated, as described below. 

In this paper we propose an algorithm to produce uncorrelated jumps along 
each path, without generating such a multitude of Markov chains. This method 
is based on stochastic step functions (SF), which will be defined below. As 
opposed to MH, it is not an accept/reject algorithm. Therefore it is able to 
generate an uncorrelated chain of values distributed according to any given 
probability measure. We will test this algorithm in jump-diffusion computer 
simulations and compare with MH. We apply the simulation techniques to a 
Gaussian process, for which the exact distribution at t = 1 is known. Conver- 
gence towards this exact result is studied. 

We also consider adaptive variants of these algorithms. For these, we sub- 
divide K into appropriate regions, and generate a discrete distribution that 
allows us to efficiently draw among these regions in an uncorrelated way. As 
the simulation progresses, this discrete distribution is adaptively improved. The 
calculation of the jump rate A is also adaptively improved. 

As an interesting application of these simulation techniques, we also look at 
infinite activity pure jump processes, which are also a Levy process subclass. 
Here, motion occurs in the form of an infinitude of discontinuous jumps. Some 
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of these processes can be approximated by jump-diffusion by substituting the 
smallest jumps for an appropriate Brownian motion [T]. The examples we focus 
on are the NIG and CGMY processes. 

In section [2] we review the elementary facts about jump-diffusion processes. 
In section [3] we describe the relevant simulation methods, and discuss their pros 
and cons, as well as provide results from numerical experiments. In section [4] we 
describe our probability space subdivision method and the accompanying ad- 
aptability properties of the algorithms, and study the efficiency and correlation 
strengths of different algorithms by simulation from a Gaussian distribution. 
We then proceed in section [5] to perform simulations of jump-diffusion pro- 
cesses. We notice how the MH correlations adversely affect the distribution of 
the simulated process, and compare the SF and MH algorithms for simulation 
of a process with a Gaussian Levy measure. Lastly, in section [6] we review the 
jump-diffusion approximation of infinite activity pure jump processes, and apply 
our simulations techniques on the infinite activity NIG and CGMY processes, 
in order to compare algorithms in these interesting cases. 



2 Jump- diffusion processes 

First we review some elementary facts about real-valued jump-diffusion Levy 
processes on a time interval [0, T]. Such a process can be expressed as a sum of 
three simple components 

Lt^ fit + Bt + Jf (1) 

The first part is a deterministic drift with rate fi, Bt is Brownian motion, and 
Jf is a compound Poisson process. The latter describes completely the discon- 
tinuities (jumps) of the paths of Lt, by means of a Levy measure v on R. We 
will several times abuse notation by writing v both for the Levy measure and 
its density. Firstly, this measure determines the jump intensity (jump rate) 

A := v{M) < 00, (2) 

of Lt- Secondly, the corresponding Levy probability measure 

vi := v/\ (3) 

determines the distribution of jump sizes. Note that we do not have E[Lt\ ~ fJ,t 
in general, although this is satisfied in our examples. 
Now, Jt can be expressed as 

Aft 

Jt^Y^^j, (4) 

i=i 

where Nt is a Poisson process of intensity A, and the jumps {Vj} are indendent 
random variables distributed according to the Levy probability measure. We 
do not simulate the Poisson process Nt directly. Instead, for each path we draw 
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the total number of jumps individually from an exponential distribution with 
average AT. Then we randomly distribute those jumps on [0, T]. We will choose 
T = 1 for our simulations. The lengths of consecutive time intervals between 
jumps will be independent and exponentially distributed, and give the correct 
jump intensity. This is described in [3]. 

Note that for a general Levy process A is in general not finite, because the 
Levy measure might diverge too strongly towards the origin. However, the 
following condition always holds. 



which restricts the strength on the divergence of v at the origin. 

The difficulty in simulation is to draw independent jump values from ■ We 
will generate a Markov chain with v\ as its invariant measure. However, it is 
essential for correct simulation that jumps along each path are uncorrelated. 
Note that existence of correlations between jumps in different paths will not 
ruin the convergence in distribution, but only slow it down. 

For a Markov chain {Xi}, we define the sequential correlation as 



where X is the Markov chain average and cP' is its standard deviation. 



3 Simulation of a probability measure 

To simulate the jumps, we must draw independent values from v\. In cases 
where is complicated, it is common to use the Metropolis/Hastings (MH) 
algorithm jTJSI- This method generates a Markov chain with as its invariant 
density. Unfortunately, the Markov chain often has large correlations between 
successive values. Successive values in such a chain cannot be used to represent 
jumps within a single jump-diffusion path. 

It is possible to reduce this problem by several methods. One is to skip 
a number of terms in the Markov chain to reduce correlation. This results 
in a loss of efficiency. Another method is to generate multiple independent 
Markov chains. Each jump-diffusion path then uses values from different Markov 
chains. This will lead to a correct pathwise simulation, and therefore correct 
convergence in distribution. The correlations will in this case only slow down 
the convergence. 

We will look at two methods, MH and one based on stochastic step fuctions 
(SF). Both rely on a transition probability distribution p to determine upcoming 
values from the previous one. We say that we are dealing with a local algorithm 
if p is localized around the current value. An independent algorithm results if p 
is independent of the current position. 




(5) 



(6) 



X :^ E[X], a^:^E[X^-X\ 
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Use of an independent MH algorithm can reduce correlations dramatically 
if one is able to find a p that approximates vi . The downside is that efficiency 
will suffer if p is not computationally simple. We will focus on generation of low 
correlation Markov chains in order to get by using only one chain for the Levy 
process jumps. 

3.1 Metropolis/Hastings (MH) 

The well-known MH algorithm generates a realisation {xi} of a Markov chain 
distributed according to vi. Its most popular incarnation is local, where the 
transition probability density p{xi^Xi+i) involved in each transition Xj^ I y Xi^\ 
is localized around the value Xi, and its width is adjusted to give an average 
acceptance rate somewhere around 1/2. The resulting correlation between suc- 
cessive values can be classified into two types: 

• Rejection correlation: Since both the local and independent algorithms are 
based on acceptance/rejectance, correlations are introduced by repeated 
values due to rejections. With an acceptance rate around 1/2, repeated 
values will often occur. 

• Transition correlation: In order for the local algorithm to be efficient, 
its transition probability measure must in many cases be quite 
strongly localized; otherwise too small an acceptance rate will result. Thus 
the difference between subsequent variates of the Markov chain will tend 
to be small. 

To reduce transition correlation, the width of the transition distribution can be 
increased. But this typically leads to a reduced acceptance ratio, which gives 
an increased rejection correlation, and vice versa. 

As a numerical example, consider a unit variance normal distribution, simu- 
lated with a simple localized uniform transition probability measure. The correl- 
ation defined in (|6| as a function of transition probability width is displayed in 
Figure^ As expected, the correlation has a minimum. Towards smaller widths 
the transition correlation increases, and towards larger widths the rejection cor- 
relation increases. Thus there is a lower bound on the amount of correlation for 
this algorithm, and it seems that the local MH algorithm is therefore not suited 
for our purposes. We will from now on focus on the independent MH algorithm. 



3.2 Stochastic step function (SF) 

Next, we propose a simulation algorithm based on stochastic step functions, 
that can be adjusted to completely eliminate correlations. This possibility of 
vanishing correlations is an attractive property that allows it to be used as a 
reference algorithm. It can also be adjusted to run more efficiently, with nonzero 
correlation. 
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Figure 1: Correlation c versus transition probability width w for the local MH al- 
gorithm, in a simulation of the unit normal distribution using a uniform centered 
transition probability distribution. Correlation is minimal around w ~ 7, where 
c w 0.55. 

Let us give a pathwise definition of the SF process on [0, cxd) for some meas- 
ure density on f2 C M, not necessarily normalized. Consider a sequence 
{tJ^q C (0, oo) of resting times and corresponding jwm^> times {tij^^O' defined 
recursively by 

to = 0, := ti + Ti- 

In addition, let {si}^o C O be a Markov chain with initial value sq = and 

transition probability density /9(.Sj, Si+i). Assume that p{si, •) is absolutely con- 
tinuous with respect to the Lebesgue measure on M, and homogeneous, i.e. it 
can be expressed as p{si+i — Si). 

From this data, we are now ready to define our stochastic step function 
process Xt, by defining its paths as the piecewise constant cadlag functions 
(right-continuous with left limits) of the form 

oo 

Xt = ^ SiXh{i), h ■= [ti,ti+i), 

where xi is the indicator function for the interval /, and where the resting times 
are given by Tj := v{si). 

We see that the paths of Xt linger for a some time at each of its attained 
positions, with a resting time equal to the local value of the density v. Con- 
sider the graph of such a path on [0,T] for large T compared to sup!/. As T 
increases, the relative path length within a given subset C on the vertical 
axis converges towards v(^A)lv{^ =: v\{^A). When sampling this path on a 
uniform time-grid, the obtained values will therefore be distributed according 
to the probability density v\. 
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As in the case of local MH, we get transition probability correlations if p 
is localized. In this algorithm, however, there is no accept/reject step, and 
therefore no rejection correlation. As an example, let us choose p to be the 
uniform probability distribution on VL. As noted above, if we sample the step 
function path generated above on a uniform time grid, we get a chain of values 
distributed according to vi . If we make sure that the time discretization interval 
size is larger than sup v, repeated values will never occur and there are no 
correlations. 

Note that as in MH algorithms, we do not need to know scale factor relating 
vi and V. If sup v is initially unknown, we can simply start with any estimate, 
and improve it adaptively as the algorithm runs. 

It is easy to see that the amount of computer time spent between each 
discrete time value is proportional to supz//z/(r2) in its correlation-free mode. If 
this ratio is large, the algorithm will be inefficient. 

One can make the time discretization finer to increase the rate of variate 
generation, but this introduces correlations since values for which v is large will 
be repeated with complete certainty. This is similar to the Metropolis algorithm, 
where values of large v are repeated (although not with complete certainty) . The 
SF algorithm can be made less deterministic in this case by letting the resting 
times be random variables distributed according to an exponential distribution 
with mean v{x), where x is the current position. 

Since we are concerned with minimizing correlations in the context of jump- 
diffusion simulations, we will use independent algorithms, where is 
independent of Xi . 



3.3 Simulation comparison 

Let us illustrate the advantages of the local SF algorithm over local MH by 
considering an example with a probability distribution with several modes on 
a sample space Vl — [0,1]. Our goal is to simulate values from a probability 
distribution proportional to the following unnormalized density with two strong 
modes, 

r 1 ,a;e [0,0.25) U [0.5, 0.75) 
^ ' \ 0.01 , otherwise. ^ ' 

For both MH and SF, we used a localized uniform transition probability density 
p of width 1/2. For MH, this gave an acceptance rate of approximatively 0.55. 
The SF Markov chain realisation was obtained by sampling the stochastic step 
function using a lattice spacing slightly larger than sup v in order to avoid any 
repeating values. 

The results are shown in Figure |2] One sees immediately that the MH 
algorithm has the potential of getting stuck inside a mode for long periods. 
This is caused by rejection correlation as found in accept / reject algorithms such 
as MH. On the other hand, the SF algorithm has no such problems since it lacks 
such correlations. This serves to illustrate a problem that often can occur with 
MH. The SF algorithm is a clear favourite in this case if mixing is important 
for the application of these Markov chains. 
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Figure 2: Markov chain realisations using local MH and local SF. The SF al- 
gorithm mixes much better than MH. 



4 Probability space subdivision and adaptability 

The Levy probability measure from which we will simulate in the context of 
jump-diffusions will often be defined on the unbounded probability space 17 = M. 
This presents no difficulty for the MH algorithm. SF algorithms on the other 
hand need to impose upper and lower cutoffs on Q, in order to avoid step 
functions that diverge towards infinity. For simplicity, we impose such cutoffs on 
both algorithms in our examples. Alternatively, one could perform a topology- 
changing coordinate transformation on Q to obtain a compact space, however 
we will not do this in our examples. From now on we therefore assume C M 
to be a bounded interval which we choose symmetrically about the origin. We 
will only deal with symmetrical Levy measures in our examples. 

In order to reduce correlations in the MH algorithms, and increase efficiency 
in the SF algorithms, we define a finite disjoint subdivision {Ui} of SI, where 
UiUi = SI. We construct a discrete distribution i> on the finite set {Ui}. For each 
Ui C S7, the value i){Ui) must approximate i^{Ui), i.e. the Levy measure weight 
of Ui. The simulation algorithm starts with a preliminary estimate of i> which is 
adaptively improved throughout the simulation, as explained below in the two 
cases of MH and SF. This is reminiscent of variance reduction schemes used in 
Monte Carlo integration, such as stratification and the VEGAS algorithm |5]. 
We now describe in more details how this is done in each case. 

4.1 Adaptive Independent MH (AIMH) 

As explained above, we will use an independent MH algorithm. The independent 
transition probability p is defined as follows. First, use the discrete unnormalized 
probability measure i> to draw a subset Ui. Then a random position within 
this subdomain is proposed and the value of i/ at this position is calculated. 
Thereafter follows the usual MH accept/reject step. 
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The initial draw of Ui is done without introducing correlations, using e.g. 
an efficient algorithm which is independent of the normalization of i/ [S]. The 
registered values of v are accumulated, and used periodically in the simulation to 
improve j>. Essentially, a separate Monte Carlo simulation is being performed 
within each subdomain Ui to improve the discrete distribution i^, while the 
algorithm proceedes to generate new draws from v. 

Since subdomains U of are drawn by means of i> without introducing cor- 
relations, the amount of correlation generated by the algorithm as a whole is re- 
duced. Transition correlation is completely eliminated since since the algorithm 
is independent. It is impossible to completely eliminate rejection correlation in 
an MH algorithm unless p is identical to vi . However, since u is well approxim- 
ated on each subdomain (as long as the subdivision is sufficiently fine), these will 
be small, and the acceptance rate will be high. As z/ adapts, the acceptance rate 
rises and correlation decreases. The subdivision discretization implies a lower 
bound for the amount of correlation. But this is still a big improvement on a 
basic local MH algorithm for the cases we consider. Note that the subdivision 
cannot be made arbitrarily fine, since the discrete algorithm will then become 
inefficient. 

4.2 Adaptive Independent SF (AISF) 

It is possible to improve the efficiency of the SF algorithm by a similar construc- 
tion, without introducing correlations. In this case, the SF algorithm proceeds 
as follows. First draw a subdomain Ui using i). Draw a random position x 
within this subdomain, and record the value v{x) for later use to improve v. 
For each subset Ui we keep an estimate of supj v which is updated for each 
sample. Each subset is also accompanied by a local time variable ti. For each 
position X within Ui, we increase ti by the resting time i'(x). We choose new 
positions independently within Ui until ti exceeds the current estimate of supj v. 

When the above loop exists, we have determined our new position and we 
subtract supj v from ti in anticipation of the next time we enter Ui. In a sense 
we are using a different time discretization in each subset Ui of il. 

This modified algorithm avoids spending a lot of time in areas of low prob- 
ability for two reasons. First, the low probability subdomains will rarely be 
entered when drawing from the dicrete distribution i). Secondly, the amount 
of time necessary to exit the loop in those regions is much smaller, since the 
local supi^ is small. By exploiting the subdivision of fl, and using an efficient 
algorithm for drawing from discrete distributions, we improve the efficiency of 
the SF algorithm drastically in nontrivial cases. 

Both i) and the estimates sup v are adaptively improved throughout the 
simulation. 

4.3 Simulation test 

To check the rate of convergence of the differerent methods, we simulated from 
a Gaussian distribution with unit variance. The rate of variate generation. 
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Figure 3: Histogram error versus computer time, for simulation from normal 
distribution. 



correlation, and the convergence towards the known exact result were analyzed. 

A histogram of the simulated values is gathered, and compared with a histo- 
gram calculated from the Gaussian distribution. We define the histogram error 
using an L°° measure, i.e. as the absolute value of the maximum deviation from 
the exact result. The plots show the histogram error as a function of simulation 
running time. The simulations were run on one core of an Intel T4400 laptop 
processor, but only the relative efficiencies matter here. Results are shown in 
Figure [3] and are discussed below. 



4.3.1 Local MH 

Since a low correlation will be important for our later use on Levy processes, we 
selected the transition probability width 0.75, giving the minimal correlation of 
approximatively 0.55, according to Figure [T] The variate generation rate was 
approximately 3.4 • 10^/s in this case. 
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4.3.2 AIMH 

As expected, as great improvement of the correlation was obtained compared to 
local MH. Since the algorithm is more complicated, the variate generation rate 
decreased to 2.3-10^/s, or roughly 60% of local MH. However, the correlation was 
reduced to around 0.03. This dramatic decrease results in a faster histogram 
convergence in terms of computer time. The lower correlation nature of this 
algorithm will be noticeably beneficial when simulating Levy processes. 

4.3.3 Local SF 

We adjusted the basic SF algorithm parameters to give zero correlation, and 
used a uniform transition probability density on [—5,5]. Thus we are running 
it quite inefficiently as regards variate generation speed, which turned out to be 
around 1.4 • 10^/s. Despite this, the histogram converges faster than local MH, 
due to the lack of correlation. 

4.3.4 AISF 

As expected, when including the domain subdivision and adaptive algorithm, 
the SF algorithm efficiency increased (in fact doubled) with a variate generation 
rate of 3.1 • 10^/s. The amount of correlation here is still zero, which leads to 
the fastest histogram convergence. So this is a great improvement at no cost 
other than increased algorithm complexity. It is perfectly suited for simulating 
consecutive jumps in jump-diffusion process paths. 

5 Jump- diffusion simulation 

As previously mentioned, existence of correlations among jumps within the same 
jump-diffusion path is disastrous. Let us check this by using an local MH 
algorithm to simulate an exactly solvable stochastic process [6] and compare 
distributions at i = 1. It is defined as 



where Nt is a Poisson process of intensity A. We choose a — 1, — 0, S ~ 1, 
and A = 10 in order to get an appreciable number of jumps within the time 
domain [0, 1]. 

The exact PDF of this process for i > is [3, Eq.(4.12)] 



i=l 




.2 
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Figure 4: Local MH simulation of distribution at t = 1 of Gaussian process, 
and exact result. Due to jump correlations, the simulated histogram has heavy 
tails. 



For local MH simulations, the positive correlations between subsequent jumps 
leads to heavy distribution tails. Results confirming this are shown in Figure |4j 
The simulation of the jumps used a localized uniform transition probability of 
width 4, which resulted in a MH acceptance rate of around 0.56. 

We now turn to our serious attemps at jump-diffusion simulation. We have 
simulated the same process using AIMH and AISF algorithms. In this case we 
use the already known value of A = 10 in the simulations, so the adaptability 
consists solely of the calculation of the discrete subdivision distribution, and in 
the AISF case also of the local subdivision supremum calculations. As opposed 
local MH, in this case the distribution at i = 1 approaches the exact curve, 
as seen in Figure [s] Note that despite our use of upper/lower cutoffs in the 
implementation of the discrete subdivision, the tails of the local MH simulation 
are still somewhat heavy. 

We use the same L°° measure of convergence detailed above. Convergence 
results are shown in Figure [6j The AISF outperforms the AIMH algorithm due 
to its lower correlation which leads to faster and more accurate convergence. 



6 Application to NIG and CGMY infinite activ- 
ity pure jump processes 

We present two applications within the realm of infinite activity pure jump Levy 
processes, namely NIG and CGMY. To this end, we employ the jump difffusion 
approximation of these processes [1]. 
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Figure 5: Distributions of Gaussian process at t = 1 for AIMH and AISF 
simulations. The continuous curve represents the exact result. Notice the slight 
heavy tails in the AIMH case, due to nonzero jump correlations. 
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Figure 6: Error in distribution at f = 1 for the Gaussian process, versus com- 
puter simulation time. AISF gives the fastest convergence. The AIMH heavy 
tails cause a lower bound on the error. 
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6.1 Jump-diffusion approximation 

For an infinite activity Levy process, the intensity A is undefined, since its defin- 
ing integral ([2| diverges. However, by virtue of (|5|, it is possible to approximate 
the infinitude of small jumps by a Brownian motion process. Consider an infin- 
ite activity pure jump Levy process Lf. For e > 0, define the following subsets 
of f7, 

A^'- := {\x\<e:xeM.} 
A'^+ := {\x\> e:x eR}. 

These represent small and large jumps, respectively. 

We now define the unique Levy measures v'^'~ and on these subdomains 
as follows. For any z^-measurable U cM. — {0}, define 

v''+{U) u{Ur\A''^+). 

Each of these Levy measures in turn define its own Levy process LI and L*^ of 
infinite and finite activity, respectively. We have the unique process decompos- 
ition 

Lt^Lr+Ll'+. (8) 
It follows from ^ that the intensity of the large jump component process, 

is well-defined, and so L^'^ is a compound Poisson process. For small e, the small 
jump process Ll'~ can be often be approximated by the following Brownian 
motion [1] 

Ll'-^a{e)Wu ^(e)^ := Var[Lr], (9) 

where Wt is a Wiener process. The variance must exist for the approximation 
to be well-defined. A sufficient condition is [T] 

lim ^ = oo. (10) 

c-i-o e 

Thus we have the approximation 

Lt^a{t)Wt + L'^+. (11) 

Any additional nonvanishing deterministic drift component of Lt would appear 
trivially on the right hand side of this equation. 

In these cases, the intensity A is not known beforehand, and also depends 
on our choice of e. Since our algorithms already calculate v, and the individual 
volumes of each subset in the subdivision of Q, is known, it is easy to use this 
to update an estimate for A concurrently. 
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6.2 Simulation of NIG 



We can now check the quaUty of the jump-diffusion approximation in conjunc- 
tion with our AIMH and AISF algorithms. Since the density of NIG is known, 
and a direct algorithm [3 , Alg.6.12] exists for simulating from its distribution at 
any time, we have ample material to compare our results. The numerical results 
for the density at t = 1 are collected using logarithmic plots in Figure [7] in order 
to emphasise the distributional tail behaviour. In terms of the parametrisation 
used in (3], we use a = 1, 9 — 0, n — 1/2. Figure [t] contains the results from 
the direct simulation algorithm, as well as the results of the jump-diffusion ap- 
proximation using AIMH and AISF, where we have used the small jump cutoff 
value e = 0.005. 

The volatility a of the Brownian component in the jump-diffusion approxim- 
ation, defined by (|9|, is calculated analytically using an approximate expression 
for the Levy measure at small The Levy jump density is 

,y{x) = ^e^-K,{a\x\), (12) 
7r|a;| 

where Ki is the irregular modified cylindrical Bessel function of first order. As 



an asymptotic approximation of (12 1 at small \x\, we use 



e^'^^l, K,ia\x\)^^. 

a\x\ 



Inserting the chosen parameter values, and using (|9| which gives ct^ as the 
second moment of the Levy jump density on [— e, e], we get 

a w V2e(5/7r w 0.067. 



Notice that this expression satisfies the Asmussen/Rosinski condition ( 10 ). Ow- 
ing to this small value, the Brownian part of the Levy process has a negligible 
influence on the distribution tails at t ~ 1, compared to the contributions from 
larger jumps. 

The direct algorithm from [3J works best, as expected. We are doing this 
to compare AIMH and AISF, for a known case with an exact solution. The 
AIMH and AISF algorithms are general and can be used on any Levy process 
that allows a jump-diffusion approximation. This case gives us an indication of 
how trustworthy these methods are when applied to more exotic Levy processes 
for which we do not have an exact result or a simplified simulation methods as 
in this case. 

Most sources of errors are common to both algorithms. These are related 
to inaccuracy in the calculation of the a and A, the latter being related to the 
cutoff imposed on the jump domain fi. This cutoff will cause an inaccuracy in A 
since the total weight of ft will not be calculated. Improved subdivision schemes 
of Q are possible, e.g. employing coordinate transformations that transform D, 
into a compact domain. We have not done this, since it unimportant for the 
matters we are considering. 
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Figure 7: Simulated and exact NIG PDFs at i = 1. The AIMH algorithm has 
heavy tails. 
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Figure 8: NIG histogram error at t = 1, versus computer time. The AISF 
clearly converges faster. AIMH reaches a steady state in the convergence plot 
due to the error caused by the nonvanishing jump correlations. 



We note that AIMH will never be completely correlation- free, and will there- 
fore tend to produce heavy tails as is obvious in Figure [7] No such problem exists 
for AISF. The cutoff on ft will naturally lead to weak tails, as seen in the plot. 
This can be remedied by enlarging the cutoff value, and/or treating large jumps 
differently. 

The measurement data for convergence analysis is plotted in Figure [8] from 
which one readily sees that the AISF algorithm converges faster and more ex- 
actly. 



6.3 Simulation of CGMY 

As a second example of a pure jump infinite activity process, we turn to CGMY 
j^. In this case we have no closed form expression for the distribution. We 
do however have the characteristic function, from which the distribution can be 
obtained by means of a numerical inverse Fourier transform. We have performed 
this calculation, and used the result as a the benchmark against which our path 
simulation algorithms are compared. 
In this case, the Levy density is 

""^""'-l Ce-«l-l/|a;|i+^ ,x<0. ^^-^^ 

The parameter space for the CGMY model is C, G, Af > and Y < 2. By ex- 
panding the Levy density in negative powers of x around the origin and keeping 
only the lowest order terms, we get 
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Figure 9: CGMY PDF at i = 1, simulated and exact results. Note again the 
heavy tails in the AIMH case. 



so by ( 10 1, the jump-diffusion approximation is valid only for < F < 1. In fact, 
the volatility can in this case be expressed exactly in terms of the incomplete 
gamma function, which we will not show here. 

We simulated the CGMY process with parameters C = G = M = l, Y — 
1/2, using e = 0.005, and 100000 paths. In this case, the volatility for our choice 
of parameter values is 

cr« 0.022. 

Also in this case it has negligible influence on the distribution tails at t = 1. 

The results for the distribution at i = 1 of the process is given in Figure [9] 
and the convergence measurements are plotted in Figure [TOj The conclusions 
are similar to the NIG case. 



7 Conclusions 

We have studied two algorithms for jump-diffusion and infinite activity pure 
jump process simulation via jump-diffusion approximations. Most importantly, 
we have studied the proposed SF algorithm based on a stochastic step function. 
It has some advantages over MH accept /reject algorithms. It is possible to 
configure it to have an arbitrarily small autocorrelation. Our simulations show 
that in the case of simulation of Levy processes, this algorithm can represent 
an improvement over the MH method that we have considered. The numerical 
results show an improvement in the tails of the distribution of the Levy process 
at a given time while at the same time converging faster. The MH algorithm 
will tend to give heavy tails, due to the problem of positive correlations between 
large jump values. 

In our computer simulations, we also implemented a subdomain discretiza- 
tion with a corresponding adaptively improved discrete probability distribution. 
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Figure 10: CGMY histogram error at f = 1, versus computer time. 



This method helps to reduce correlations for the MH algorithms, since the sub- 
domains themselves are drawn without using an accept/reject algorithm. In 
the SF case, it improves the variate generation speed while maintaining the 
correlation-free nature of the method. 
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